www.gusucode.com > muon generator源码程序 > muon generator源码程序/Geant4_Muon_histogram_generation_v1.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Geant4-MATLAB muon generator file (Version 1).
%
% This MATLAB script has been developed to facilitate the generation of the
% muon energy and angular histograms use with the Geant4 macro file.
% The script generates a look-up table that contains the sampled muon
% angles and energies as well as histograms for muon angular and energy distributions.
%
% This muon generation is based on a methodology that combines the
% Smith & Duller (1959) phenomenological model and statistical sampling
% algorithms. The inverse transform provides a means to generate samples from the
% muon angular distribution. The Acceptance-Rejection method is used to
% provide the energy component. Final output provides the user defined
% histograms for use in the Geant4 macro file.
%
% For non-standard energy and angular distributions Geant4 requires the
% utilization of the General Particle Source module which via the
% G4GeneralParticleSource class allows specification of user defined
% angular and energy distributions. The user defined histograms are
% specified in macro files using the following commands:
%
% 1. /gps/particle mu+ (specify particle type)
% 2. /gps/ang/type user (user defined histogram)
% 3. /gps/hist/type theta (zenith angle histogram)
% 4. /gps/hist/point Bt Wt (angular histogram values)
% 5. /gps/ene/type Arb (user defined histogram)
% 6. /gps/hist/type arb (point-wise energy spectrum)
% 7. /gps/hist/point Eh Hh (energy spectrum values)
% 8. /gps/hist/inter Lin (interpolation scheme: Linear)
% 9. /run/beamOn 1000 (number of particles)
%
% where a short explanation of each command appears in parentheses.
% The histograms represent differential functions and must be included
% one bin at a time. Angular histogram is described using the bin upper
% boundary and the area of the bin. Energy spectrum (point-wise) is
% described using the bin center and the height of the bin. The first
% value of each histogram must be the lower boundary of the bin and a
% dummy value that is not used.
%
% Instructions: Run the file. On the command line set the minimum and
% maximum muon energy and the minimum and maximum zenith angle.
% Output:
%
% 1. The output "Muon_table" matrix contains the angles and energies
% of the sampled muons
% 2. The "User_defined_hist_energy" matrix contains the point-wise energy
% (MeV) spectrum for the Geant4 macro file. Just copy paste to your Geant4
% macro file.
% 3. The "User_defined_hist_angle" matrix contains the zenith
% angle (radians) histogram for the Geant4 macro file. Just copy paste to your
% Geant4 macro file.
%
% This file is free for use. More details, examples and validation results
% can be found on the journal papers:
%
% S. Chatzidakis, S. Chrysikopoulou, L.H. Tsoukalas (2015)
% "Developing a cosmic ray muon sampling capability for muon tomography
% and monitoring applications", Annals of Nuclear Energy, to appear
%
% S. Chatzidakis, L.H. Tsoukalas (2015)
% "A Geant4-MATLAB muon generator for Monte-Carlo simulations",
% Transactions of American Nuclear Society, Vol. 113, to appear
%
% Users are kindly requested to cite the above journal papers in their
% publications. Please comment. Thank you!
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Set minimum and maximum muon energy (Validated range 1-100 GeV)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Emin=input('Enter minimum muon energy (GeV):');
Emax=input('Enter maximum muon energy (GeV):');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Set minimum and maximum muon zenith angle (Range 0-90o)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Theta_min=input('Enter minimum muon zenith angle (degrees):');
Theta_max=input('Enter maximum muon zenith angle (degrees):');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Set the desirable number of muons to be sampled
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N=100000;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Phenomenological model constants
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Alpha=0.002382; % constant A
lambda=120; % absorption mean free path 120 g/cm2
kappa=2.645; % exponent (-)
bp=0.771; jp=148.16; % correction factor (-); factor (GeV)
alpha = 0.0025; % muon energy loss in GeV/g/cm2
rho = 0.76; % fraction of pion energy that is transferred to muon
y0=1000; % atmoshperic depth g/cm2
bm=0.8;Bm=1.041231831; % correction factor (-);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Inverse trasform and Accept-Reject method
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
q=1; % initialize
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Inverse trasform to sample muon zenith angle
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
theta1=Theta_min*pi/180;theta2=Theta_max*pi/180;
for i=1:N
tm=theta1+(theta2-theta1).*rand();
tm1=tm/(pi/2);
theta(i)=acos((1-tm1)^(1/3)); % select muon angle (in radians) using inverse transform
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Calculate maximum value of phenomenological model
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Em=[Emin:0.1:Emax];
Ep=(Em+alpha*y0*(sec(theta(i))-0.1))*(1/rho);
Pm=(0.1*cos(theta(i)).*(1-(alpha*(y0*sec(theta(i))-100)./(rho.*Ep)))).^(Bm./((rho.*Ep+100*alpha)*cos(theta(i))));
f=Alpha.*Pm.*(Ep.^(-kappa))*lambda*bp*jp./(Ep.*cos(theta(i))+bp*jp);
C3=trapz(Em,f); % normalization constant for f
C4=max(f); % maximum value of phenomenological model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Accept-Reject method to sample muon energy
% (for the previously selected muon angle theta)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k=1:1000000
Em1=(Emin+(Emax-Emin)*rand()); % select a random energy from uniform U(Emin,Emax)
u=rand(); % select a random number from uniform U(0,1)
Ep1=(Em1+alpha*y0*(sec(theta(i))-0.1))*(1/rho); % calculate pion energy
Pm1=(0.1*cos(theta(i))*(1-(alpha*(y0*sec(theta(i))-100)/(rho*Ep1))))^(Bm/((rho*Ep1+100*alpha)*cos(theta(i)))); % calculate probability
f1=(1/C3)*Alpha*Pm1*(Ep1^(-kappa))*lambda*bp*jp/(Ep1*cos(theta(i))+bp*jp); % calculate intensity based on sampled energy and angle
f2=C4/C3;
f3=f1/f2; % f(y)/g(y)=f1/f2 where g(y)>f(y)
if u<=f3 % if u<=f(y)/g(y) then set X=Y
Muon_Energy(q)=Em1; % accept energy
q=q+1;
break
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Create a lookup table: 1st column is sampled muon angle (radians),
% second column is sampled muon energy (GeV)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Muon_table=[theta;Muon_Energy]';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Set bin size and produce histograms for plotting
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
bin1=2*(numel(Muon_Energy)^(1/3)); % bin size for energies selected according to Rice rule
[C,x2]=hist(Muon_Energy,round(bin1));
C1=C/trapz(x2,C); % normalization of the histogram
bin8=2*(numel(theta)^(1/3)); % bin size for angles selected according to Rice rule
[C8,x8]=hist(theta,round(bin8));
C9=C8/trapz(x8,C8); % normalization of the histogram
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Calculate bin size and bin area for angular histogram
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
bin_size=(max(theta)-min(theta))/numel(x8); % bin size
C10(1)=0; % dummy value
C10(2:(numel(x8)+1))=bin_size.*C9; % area of bin
x9(1)=min(theta); % first bin lower boundary
x9(2:(numel(x8)+1))=x8+(bin_size/2); % bin upper boundary
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Calculate bin center and height for point-wise energy spectrum
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x3(1)=x2(1); % first bin value
x3(2:(numel(x2)+1))=x2; % bin centers
C2(1)=0; % dummy value
C2(2:(numel(x2)+1))=C1; % bin height
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Generate commands for use in Geant4 macro files
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C3=repmat(['/gps/hist/point'], numel(x3),1);
C11=repmat(['/gps/hist/point'], numel(x9),1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Generate energy point wise spectrum (in MeV)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
User_defined_hist_energy=[char(C3) blanks(numel(x3'))' num2str(x3'.*1000) blanks(numel(C2'))' num2str(C2')];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Generate zenith angle histogram (in radians)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
User_defined_hist_angle=[char(C11) blanks(numel(x9'))' num2str(x9') blanks(numel(C10'))' num2str(C10')];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Plot histograms
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
plot(x8,C9,'o') % plot distribution of zenith angles
xlabel('Zenith angle (rad)')
ylabel('Arbitrary units')
figure(2)
loglog(x2,C1,'o') % plot distribution of energies
xlabel('Energy (GeV)')
ylabel('Arbitrary units')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Clear workspace
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clearvars -except Muon_table User_defined_hist_energy User_defined_hist_angle % clear workspace